function conditionNumber = calculateConditionNumber( preconditioner, ErrTol, MaxIts, verbose)
%CALCULATECONDITIONNUMBER Summary of this function goes here
%   Detailed explanation goes here
    n=size(preconditioner,2)-1;
    H=makeHFilter2D(n,preconditioner);

    j=10; %hardcoded... yak
    W=1./(2*(rand(n)-0.5)); % W is distributed N(0,1) iid
    W=W-sign(W);
    W=W/norm(W(:));% and normalized
          L=3*n-2;

        Hhat=fftn(H.filter,[L L]);  
    if 0

        temp=W;
        for i=1:j
            temp=applyHFilter2D(temp,Hhat);
            if i == j-1
                last=temp;
            end
        end
    end
    v=W;
    rho=1;
    for j=1:j
        v=feval(@applyHFilter2D,v,Hhat);
        d=norm(v(:));
        rho=sqrt(d);
        v=v/d;
        %fprintf('iter=%d\trho=%e\n',j,rho);
    end
    
    %normH=sqrt(norm(temp(:))/norm(last(:))); %upto factor 10

    normH=rho;
    
    if nargin<4
       verbose = 0;
    end

    if nargin<3
       MaxIts = 10;
    end

    if nargin<2
       ErrTol = 1.e-2;
    end
    
    %temp=W;
    
    if 0    
        for i=1:j
            [temp,flag,residual,iter] = ConjugateGradient('applyHFilter2D',temp,{Hhat},ErrTol,MaxIts,zeros(size(temp)),verbose);
            if i == j-1
                last=temp;
            end
            if flag
               fprintf('Inversion did not converge. Residual error %-2.5e, iter = %d\n',residual,iter);
            end
        end
    end
    %normInverseH=sqrt(norm(temp(:))/norm(last(:))); %upto factor 10

    
    v=W;
    rho=1;
    for j=1:j
        %v=feval(applyHFilter2D,{v,H});
        [v,flag,residual,iter] = ConjugateGradient('applyHFilter2D',v,{Hhat},ErrTol,MaxIts,zeros(size(v)),verbose);
        d=norm(v(:));
        rho=sqrt(d);
        v=v/d;
        if flag
           %fprintf('Inversion did not converge. Residual error %-2.5e, iter = %d\n',residual,iter);
           %fprintf('iter=%d\trho=%e\n',j,rho);
        end
    end
    normInverseH=rho;
    conditionNumber=normH*normInverseH;

end

